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A numerical algorithm based on the finite element 
method of analysis of the boundary value problem in a 
continuum is presented, in the case where the plastic 
response of the material is given in the context of en- 
dochronic plasticity* The relevant constitutive equation 
is expressed in incremental form and plastic effects are 
accounted for by the method of an induced pseudo-force 
in the matrix equations. 

The results of the analysis are compared with observed 
in the case of a plate with two symmetric notches 
and loaded longitudinally in its own plane. The agreement 
between theory and experiment is excellent. 

INTRODUCTION 

The greatest difficulty encountered in the application 
of the classical theory of plasticity is the lack of 
knowledge of the configuration of the subsequent yield 
surface for the particular material at hand, and the 
experimental difficulties encountered in finding it in the 
fully three dimensional case. More importantly, however. 
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it has been observed by many experimenters that the shape 
of the subsequent yield surface and its position in stress 
space depends very strongly on the definition of the yield 
point, particularly in situations following prior deform- 
ation [1-3] . 

The essential premise of the classical plasticity 
theory is the assumption of an a priori existence of a 
yield surface. This implies a finite elastic domain. 

From the mathematical standpoint, a finite domain is 
necessary because of the requirement that the increment 
in plastic strain be normal to the yield surface. Thus, 
the direction of the plastic strain increment is dictated 
by the yield surface configuration. 

If plastic effects were to begin immediately upon 
loading, perforce, the domain of the yield surface would 
collapse to a point, thus making the direction of the 
plastic strain increment indeterminate since all directions 
are normal to a point. Thus, the classical plasticity 
theory cannot deal with materials that yield immediately 
upon loading. There are other difficulties associated with 
experimental attempts to describe and analyze a two-or 
three-dimensional response of a material [4] . For instance 
investigations in the hardening rule are much discussed in 
the current literature, but definitive functional forms out 
side the Prager- Ziegler rule are very few, and lack firm 
experimental verification. This rule specifically can have 
only limited application, and is inappropriate for 
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complicated loading histories. Moreover, it gives rise to 
large discrepancies between calculated and experimental data 
in loading-unloading processes [1] . Other numerical 
^^icuj.'fcies arise from the fact that the loading increments 
cannot be assigned arbitrarily a priori. When the current 
loading increment makes the stress state of a particular element 
traverse the yield surface it is necessary to come back to 
the preyious loading state and adjust the magnitude of the 
new increment of loading to ensure that the new stress state 
is located just on the yield surface. Certainly, this process 
increases the time of computation. 

In 1971, Valanis proposed an alternative theory of 
viscoplasticity called "endochronic theory" [5,6], which 
is based on irreversible thermodynamics and the concept of 
intrinsic time. The theory provides a unified point of view 
to describe the elastic— plastic behavior of materials since 
it places no requirement for a yield surface and a "loading 
function" to distinguish between loading and unloading. 

In a series of recent works, Valanis, Wu and others 
demonstrated that the endochronic theory could apply 
®ore precisely to situations involving unloading and cyclic 

behavior of metals, as well as wave propagation in the plastic 
region. 

However , in all of the works , involving more than one 
dimension, where the loading was quasi-static, the stress 
fields were homogeneous. 


129 



In the present paper a numerical algorithm is first 
implemented in a computer program, which can be used to 
analyze the material response in monotonic and cyclic loading 
in the case of plane stress or plane strain* The calculated 
results are then compared with the data obtained from a 
specially designed experiment on a notched plate cyclically 
loaded in its own plane. The validity of the endocttronic 
analysis, using this numerical algorithm, is thereby 
demonstrated. 

AN INCREMENTAL FORM OF THE ENDOCHRONIC ELASTOPLASTIC 
CONSTITUTIVE EQUATION IN TERMS OF {do} and {dc} 

The following are the formulae concerning the endoch- 
ronic constitutive equations for plastically incompressible 
isotropic materials and small deformation [7] 

2 3e P 

s * / p(z-z') dz* (2.1) 

o 

d; - || de^ || (2-la) 


dz * 


dc 

fTsT 


(2.1b) 


where p(z) and f(c) are two material functions namely the kernel 
function and hardening function respectively. 
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(2.3) 


de tf * de - -s— ds 


By definition 


de. . * de. . - i de $ . . 
J-3 13 3 aa ij 


ds. . =* da. . - ■=■ da 5. . 
13 13 3 aa 13 


(2.4a) 

(2.4b) 


In this paper the form of p(z) given by equation (2.5) 
was used in equation (2.1) 


00 

p(z) = l 

-a z 
c r e . 

(2.5) 

r=l 

with the conditions that 6 

r and R r are positive 

for all 

r and 

00 

00 

n. 


l <= r * 00 

“W 

A 

8 

• 

(2. 6a, b) 

r~l 

r=*l r 



This form of p(z) is continuous and differentiable in (0,«) 
and therefore the incremental form of equation (2.1) specified 
below can be used in conjunction with a finite element code. 

Specifically in the case where the infinitely large value 
of p (0) is approximated by a suitably large value, as is done 
in this paper, one may differentiate equation (2.1) with 
respect to z to obtain the following differential form of 
the endochronic constitutive equation: 


ds » p(0)de^ + h(z)dz (2.7) 

where 

2 . 9eP 

h(z) * / p(z-z') jp- dz' (2.8) 
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and 


(2.8a) 


* if 


The elastoplastic constitutive equations (2.3) and (2.7) , 
can then be combined and expressed in the differential form 


ds ij = 2 ^ de ij + h i; .(z)dz> 

where 

tf - p(0) {1 + 

Alternately, for computational purposes the incremental 
form given by equation (2.10) may be used, i.e., 

As ij =* 2y (Ae ij + ^oT h ij( 2 ) (2.10) 

Substituting (2.4a,b) into (2.9) and using (2.2) one obtains 
the operational incremental form of the elastoplastic con- 
stitutive equation in matrix notation as follows: 

( 2 . 11 ) 
( 2 . 12 ) 


{do} * {D> {de} + {dH p } 


where 


{D} 


C 1 c 2 


<f°2 c x 


0 
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(2.9) 


(2.9a) 


132 



and 


{dH p > 


dH 

dH 


px 


py 


dH 


pxy( 


In plane stress 


„ _ 12Ku + 4y 

C 1 *" 

x 3K + 4U 


c 2 


A A A 

6KU ~ 4y* 
3K + 4u 


« _ 2y(3k-2y) 

D, » *“ 

x 3k + 4u 


dH px = {2yh x (z) - D x h z (z) } dz/p(0) 


dH py - {2yh y (z)* - D^U)} dz/p (0) 


<*H pxy = 2yh xy (z)dz/p(0) 


In plane strain 

. 3K_+,4u 
C 1 3 


.. _ 3K - 2y 

C 2 3 


(2.13) 

(2.14) 

(2.15) 

(2.16) 

(2.17) 

(2.18) 

(2.19) 

( 2 . 20 ) 
( 2 . 21 ) 
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( 2 . 22 ) 


dH px * 2u h x (2)dz/p(0) 

a 

dH py a 2]i hy ( 2 ) dz/p (0) (2.23) 

A 

^xy = 2u ( 2 ) dz/ p ( 0 ) (2.24) 


We note that {D} is an adequate approximation to the elastic 
matrix {E}. It is evident from equation (2.9a) that when 
P ( 0 ) -*• ®, {D} becomes the elastic matrix {E}. Take plane 
stress as an example on the simple tension curve (Fig. 1) draw 


lim {D} 
p(0)-*» 


{E> = 2 (l+v) (l-v) 1 


2 2v 
2v 2 
0 0 


0 

0 

(1-v) 


■(2.25) 


We use axial tension to show the geometric meaning of equation 
(2.11). From a point A on Simple tension curve (Fig. 1) draw 
a straight line AB, the slope of which is Young's modulus 
E and its horizontal projection is de. For simple tension 


{D} {de} » Ede 


and 


BD * Ede 


(2.26) 


so BD can be considered as the first term of right hand in 
(2.11). Since CD is equal to da, the geometric meaning of 
dH is represented by the segment BC the value of which is 

Er 

negative for simple tension. 
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A FINITE ELEMENT CODE FOR THE ENDOCHRONIC THEORY OF PLASTICITY 


Using (2.11) and the principle of virtual work [11] , one 
may formulate an initial stress finite element computational 
algorithm of the endochronic theory. In fact/ we have 

fff{c) T (Se}dv * {p e3e } T ${q} (3.1) 

v ex 

and * p ex* and are respectively the vectors of nodal external 
forces and displacements of the element. Substituting (2.11) into 
(3.1) one finds that 


{K} {Aq} =• Up ex > + Upp} (3.2) 

where {K> is the stiffness matrix of the element and is the same 
as the stiffness matrix of an element in the usual elastic analysis 
but the constants C^, Cj are obtained from equations (2.14 - 2.16) 
or (2.20 - 2.21) . 

The quantity {AP p } is the incremental plastic pseudo-force 
vector for a typical triangular element used in the analysis and 
has the form 


( AP ) =* - %■ (a .AH + 8. AH ) 

< 2 i px p i pxy ; 


“Vi ' - ! (B i 4H py + a i 4H pxy ) 


i = 1/2/3 (3.3) 
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Where the components of {AHp} are given in equations (2.17 -19) 
or (2.22-24) by changing operator "d" to "A". and 8^ are 
related to the differences of nodal coordinates, i,e,, 


a i * I e ijk Ay jk ' 0 i * ’ I e ijk Ax jk (i * 1 ' 2 ' 3) 

(3.4) 

where Ay^ * y j — t Ax^ * 3C j“ x k e ijk is the P eriautation 
symbol. 

Prom equations (3.2) and (3.3) one obtains the total 
stiffness matrix {k}, total plastic pseudo force matrix 
J{AP p ) and the linear simultaneous equations for the structure. 


THE CALCULATION OF h(z) 

Equations (2.17) through (2.19) show that h(z) plays 
a central role in the calculation of (AHp) and plastic pseudo- 
force {AP >. To calculate h(z) numerically, we divide the 
P 

domain of integration (0,z) in equation (2.8) into n subregions 
whereupon 

z l^ 8e P Z i , 8e P 

h(2 m ) - / C(z m -S’) ' + / H T d*’ 

° z i-l 





p(z -z') 
in 


z 


m-1 


3e* 


(4.1) 


where z . . z . are the intial and end values , respectively , 

1— If 1 

of the intrinsic time scale of ith interval, which corresponds 
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to the ith incremental loading process, and z ffl is the current 

value of the intrinsic time scale. 

The mean value theorem, and the smoothness of e p allows 
the approximation 


- 2 n 1®£ d*' = 


I * <*»-«'» ffr <**' W 


2 * 2 . 


/ ptz^z'Jdz' 


(4.2) 

provided that there is no strain reversal in the interval 
considered. In the present work we approximate the series 
on the right hand side of equation (2.5) by three terms, i.e., 

3 -V 

P ( z) - y C e r (4.3) 


r-1 

Substituting equation (4.3) into equation (4.1) and 
using equation (4.2) we obtain the result 


h(z) 

m 


3 m 

l l C 
r=l i*l 


3eP 
r TT 


-a(z-z. .) 
e r m i-l (1 _ e 




Z»Zj 


(4.4) 


This form of h is unsuitable for numerical compulation. 

The term a r ( z ra “ 2 i«.^) ma y in the course of calculation become 
very large of ther order of 5 x 10 4 . Consequently, the 
value of the function exp{-a r ( z m “ z ^„^) } becomes a very 
small number leading to serious truncation errors. To avoid 
this difficulty we proceed as follows. By mathematical 
induction the following formula can be shown: 
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-a_Az. 


M*i> - I K« w )e, ' 

r»l 


3 

r*l 


-Gt AZi 

(e r 1 -X) 


3e p 

3Z 


z»z. 


( i 1, . .m 


(4.5) 

where h(0) ** 0 and Az. ■ z.-z. 

This is an important result to the effect that the 

history dependence of the material response (through h(z. )) at 

the intrinsic time z^ will be determined by h(z. .) and the 

x ^ x-x 


3e p 

new incremental step (through 

o Z 


and 


z=z. 


*!>• 


This formula 


is also of value in the computer program, because (a) one need 
only store the information at z^_^ to obtain results at z^, and 
(b) when using (4.5) instead of equation (4.4), the value 
of the term exp(a r Az^) is no longer small thus avoiding 
truncation errors present in the previous formulation • 

(equation 4.4) . 


THE ITERATIVE PROCESS 

For every increment of loading or unloading an initial 

value Az° is assigned to the increment of intrinsic time. 

The linear simultaneous equations are then solved and the 

displacement increments are obtained, from which the total 

deviatoric strain Ae is calculated. Also As and Ae p are 

«*• 

9 

calculated using equations (2.10) and (2.3) respectively. 

3eP 

Upon use of equations (2.1a), (2.1b) and (4.5) Az, and 

o Z 

h are obtained. Also, from equations (2.17) - (2.19) or (2.22)- 
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(2.24) AH , AH , AH and finally (AP ) are obtained. 

P* py P**Y P 

Substituting {AP^} into the simultaneous equations (4.2) we 
then obtain a new solution for the displacement increments 
as well as the other variables, including Az. The iteration 
process is continued until the difference in two consecutive 
values of Az, corresponding to two consecutive iterations, is 
less than some defined tolerance. Results are stored for 
the next step. The new loading process is then repeated. 

In this initial stress method of classical plasticity 
one [12] usually stops the iteration process if the difference 
in the magnitudes of the plastic pseudo-force vector corresponding 
to two consecutive iterations is sufficiently small. We use 
the scale Az as a criterion of convergence instead of the pseudo- 
force vector, not only because of its simplicity but because 
of its crucial role in endochronic plasticity. 

CONVERGENCE AND TOLERANCE 

The rate of convergence is very important because it 
relates to consumption of computer time, truncation error 
and other related considerations. The key of accelerating the 
convergence rate is how to choose the initial Az in order to 
begin the iteration process of a new incremental loading 
(unloading) step. An accelerator K was used to determine 
the starting value of the increment of intrinsic time Az^ by 
the relation 

Az° = K a Az i-i ’(6.1) 
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where the subscript I denotes the current incremental loading 
step and 1-1 denotes the preceding step. The superscript o 
denotes the initial value, L denotes the last value and 
is called the accelerator for the I'th increment. Equation 
(6.1) is not suitable for reversal points, at which Az° is 
taken equal to zero, because at the onset of unloading the 
response is elastic. The value of the accelerator was determined 
by the ratio of the final value of Az in the two previous 
steps, i.e., 





( 6 . 2 ) 


With the exception of the first few (three) increments 

the value of was substantially constant. To illustrate 
61 

its utility and average value of 1.24 was used and the number 

of iterations needed for convergence was compared in cases 

where K 1 = 1 and K 1 =» 0. See Pig. 2 where n pertains to 
^ sl 

the fifteenth increment and is the plastic strain near the 
tip of the notch. Curve 1 (K = 0) shosws that the convergent 

cL 

process is very slow. The reason is that at the first iteration 
Az° = 0 since K = 0 and therefore {AP}^ = 0, i.e., the loading 
process so initiated is elastic and is far away from the real 
case. Curve 2(K_ = 1) shows the convergent rate is much better 
than in curve 1, because it takes the final value of Az in the 
previous incremental loading step as the initial value of Az 
in the current step. However, in this procedure the plastic 
pseudo-load is underestimated. A value of K greater than 

cL 
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unity does increase the rate of convergence as shown in curve 
3 (K » 1.24) . Figure 3 shows the effect of accelerator 
factor K on the average iteration number N per incremental 
loading step . 

By definition the relative error ERR is defined as 

Az -Az , 

ERR = — - - n — = (4.41) 

where n is the number of iteration steps. Tolerance 
is defined as the maximum acceptable value of ERR. 

Figures 4 and 5 show the effect of tolerance on the 
accuracy and rate of convergence. In the example shown 
the smaller the tolerance the higher the accuracy (Fig. 4) , 
but the number of iterations increases (Fig. 5) . One however 
must guard against an excessively small tolerance, which may 
lie outside the inherent accuracy of the numerical computation 
and computer capability, leading to accumulation of truncation 
errors. In the present work the tolerance was 1%. 

COMPARISON BETWEEN EXPERIMENTAL DATA AND CALCULATED RESULTS 

To verify the validity of the endochronic analysis, using 
the present numerical algorithm, the distribution of strain of 
a notched specimen (made of OFHC copper) cyclicly loaded in 
its own plane was calculated and measured. One quarter of 
the specimen is shown in Fig. 6. The material functions p(z) 
and f(?) were determined by means of an experiment on a round 
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specimen of precisely the same material as the notched specimen, 
in terms of purity, grain size and treatment. The method of 
determination of these functions will not be given here but 
may be found in Ref. 13. Suffice it to say that they are of 
the following forms 

_ -a z 

P(z) - l A r e r (GPA) 

r*l 


where A^ ^ 3 * (59 2, 220, 46) and a. « _ 
x 10 ^ and 

f(?) » 1 + 0.53S 0, 72 


(27.5, 11.5, 7.67) 


The calculations were conducted on an electronic computer 
(AMDAHL 470 V/7A, close to IBM 370) in the computer center of 
the University of Cincinnati. There are 413 elements and 230 
nodes in one quarter of the specimen (Fig. 6 ) . The side of 
the smallest element is 0.25 mm. By "varying band storage" the 
amount of storage for the total stiffness matrix is 17698. 

The incremental loading for each step is 4% of the maximum 
load. The average number of iterations for each incremental 
loading was about 10 , varying from 3 to 20 . The computer time 
for each iteration was about 3.36 sec., most of which is used 
to solve the 460 simultaneous equations. The experiments were 
conducted in Metcut Research Associates Corporation. The 
strain distribution was measured using strain gauges, the 
snallest nominal length of which was 0.2 mm. Since the locations 
of the elements and the strain gauges did not coincide exactly, 
we compared the calculated results with experimental data in 
terms of plotted curves . Comparisons were made over a wide 
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range of magnitude of applied maximum stress, location and type 
of histories. 

Measured and Calculated Strain Distributions e Along the Notch 
Center Line oo' are Shown, for Applied Stress Amplitude 3.7 x 10 2 * * * * 7 PA 


(i) at first tensile peak A. Pig. 7 

(ii) at first unloading point C. Fig. 8 

(iii) at first compressive loading peak B. Fig. 9 

Letter designations as shown in those Figures. 

Measured and Calculated Strain Distributions s y Along the 

Vertical Line ob are Shown for Applied Stress Amplitude 2.3 x 10 7 PA 


(i) at first tensile peak E. Fig. 10 

(ii) at first compressive peak L. Fig. 11 

(iii) at second, loading peak H. Fig. 12 
Letter designations as shown in above figures . 

Despite the complexity of the boundary value problem and 
the inherent experimental difficulties the agreement between 
sxperimental and calculated results is excellent both from 
the aspect of tendency and magnitude. 
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Figure 1. - The Geometric Heaning of The Incremental Endochronic 
Elastoplastic Constitutive Equation 



Figure 2. - The Effect of Acceleration Factor K a on the 
Convergence of c 

^max 
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1.2 

K. 

celerabion Factor K a on the Average 
s per Incremental Loading Step 



146 


(X10 7 PA) 




Figure 6. 




(b) Fine Mesh Near the Notch Tip 
- Plate Geometry and Grid Arrangement (one quarter) 
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Figure 7. - Comparison Between Experimental and Calculated 
Distribution of Strain e y Along Notch Line 00* 

at Positive Peak A 

ft 



Figure 8. - Comparison Between Experimental and Calculated 

Residual Strain Distribution c y Along Notch Line 
00* at Point C 
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Strain Distribution c y Along Notch Lino 00' at 
Negative Peak B 



Figure 10. - Comparison Between Experimental and Calculated 
Strain Distribution e y at the Positive Peak E 

Along the Vertical Line ob at the Top of Notch 
Tip 


150 



Figure 11. - Comparison Between Experimental and Calculated 
Strain Distribution Cy at the Negative Peak L 
Along the Vertical Line ob at the Top of Notch 
Tip 



Figure 12. - Comparison Between Experimental and Calculated 
Strain Distribution ty at the Positive Peak B 
Along the Vertical Line ob at the Top Notch Tip 
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